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Abstract 



An expression is derived for the strain energy of a polymer chain under an arbitrary three- 
dimensional deformation with finite strains. For a Gaussian chain, this expression is reduced 
to the conventional Moony-Rivlin constitutive law, while for non-Gaussian chains it implies 
novel constitutive relations. Based on the three-chain approximation, explicit formulas are 
developed for the strain energy of a chain modeled as a self-avoiding random walk. In the case 
of self-avoiding chains with stretched-exponential distribution function of end-to-end vectors, 
the strain energy density of a network is described by the Ogden law with only two material 
constants. For the des Cloizeaux distribution function, the constitutive equation involves 
three adjustable parameters. The governing equations are verified by fitting observations on 
uniaxial tension, uniaxial compression and biaxial tension of elastomers. Good agreement 
is demonstrated between the experimental data and the results of numerical analysis. An 
analytical formula is derived for the ratio of the Young's modulus of a self-avoiding chain 
to that of a Gaussian chain. It is found that the elastic modulus per chain in the Ogden 
network exceeds that in a Gaussian network by a factor of three, whereas the elastic modulus 
of a chain with the generalized stretched exponential distribution function equals about half 
of the modulus of a Gaussian chain. 
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1 Introduction 



This paper is concerned with the so-called "physically motivated" constitutive equations for the 
elastic response of elastomers at finite strains. Although this subject has attracted substantial 
attention in the past half a century, it remains a focus of attention both in the communities of 
mechanical engineers, [HElISlIllElElIZlElini EDI and polymer physicists pU El EH EU H3 EH 
El E3 EH HIH US • This interest may be explained by two reasons: (i) from the standpoint of 
applications, it is induced by the necessity to have a reliable tool for the analysis of the mechanical 
behavior of polymer materials, whose response is inadequately described by conventional stress- 
strain equations, and (ii) from the point of view of fundamental research, it is driven by the lack 
of constitutive laws based on the solid ground of statistical physics that provide rather simple 
formulas for the strain energy of rubber-like materials. 

Surprisingly, the concept of Gaussian chains in an affine network introduced about a century 
ago may serve as the only example of a physically-based theory that can be employed in the 
analysis of stresses in elastomer structures (in this assertion, we exclude from the consideration 

(i) slip-link models |23| that deal with chains in non-affine networks, (ii) numerous variants, 
see [23] the the references therein, of the Kratky-Porod 1 model |25j . where the stress-strain 
relations contain the inverse Langevin function and can be resolved only numerically, and (iii) 
the Gent [22j, see also |28| . and FENE dumbbell [22] constitutive equations, which correctly 
mimic a limiting chain extensibility, but can be deduced from an appropriate statistical theory 
within the Peterlin approximation only). The model of a Gaussian chain is, however, overly 
simplified for engineering applications, because it implies the neo-Hookean expression for the 
strain energy per chain 

W = h B T(J 1 -3), (1) 

which rather poorly describes observations. Here k-Q is Boltzmann's constant, T is the absolute 
temperature, and J m stands for the mth principal invariant of the right Cauchy-Green defor- 
mation tensor C. From the physical standpoint, other shortcomings of the model of a freely 
jointed chain are that it (i) disregards short- and long-range interactions between segments, and 

(ii) implies that the end-to-end distance exceeds the contour length of a chain with a non-zero 
probability. 

To avoid these disadvantages, more sophisticated models of non-Gaussian linear chains are 
introduced. The starting point in their derivation is the treatment of a chain as a random curve 
in a three-dimensional space [3U] . The randomness of the curve reflects random fluctuations of 
the chain driven by thermal excitations. Each realization of the curve is described by the vector 
equation Q(s), where s is the arc-length. For definiteness, we suppose that s £ [0,L], where L 
stands for the contour length, and assume the end s = to be fixed at the origin, and the end 
s = L to be free. In these notation, the end-to-end vector of a chain reads Q = Q(L). 

As the description of random curves is extremely complicated from the mathematical stand- 
point, it is convenient to replace any curve by a set of N rigid segments with length I connected 
with each other. These segments are thought of as random vectors, which are independent for 
Gaussian chains or mutually dependent for more advanced models. The segment length I and 
the number of segments N are connected by the formula L = Nl. Conventionally, N is treated 
as a given integer number. The parameter I may be fixed (for the Gaussian and Kratky-Porod 
chains), or it may be considered as a random variable (for the Levy flights (SB)- For a fi xe d I, 

We employ the terminology widely used in statistical mechanics of macromolecules. In the studies on mechan- 
ics of polymers, this model is conventionally referred to as the James-Guth theory [2rJ| or the theory of Langevin 
chains |1U|. 
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the discrete model of a chain is equivalent to a random walk with iV steps that starts at the 
origin. 

Formally, a random walk is entirely described by the probability density /(Qi, . . . , Qjv)> 
where the vector Q n characterizes the nth step of the walk. For a large N, this description 
becomes, however, inconvenient, and only the distribution po(Q) °f the corresponding end-to- 
end vector Q is studied. For a spherically symmetric walk (no preferable direction in space), 
the average value of Q vanishes, Q = 0, whereas the mean-square end-to-end distance b 2 = Q 2 
characterizes the spatial dimension of a chain. The scaling analysis implies that 



b = IN", 



(2) 



where v is a scaling exponent. For a Gaussian chain without interaction between segments, 
v = ^, but for a chain that is sensitive to such interactions, v becomes higher. 

All theories of polymer chains may be divided into two groups: (i) phantom chains, for 
which two segments of a chain may be located at the same point, and (ii) self-avoiding chains, 
for which excluded volume interactions are taken into account. The first group of models is 
relatively simple, and a number of mathematically strong results have been derived for phantom 
chains. The most famous models of this group are the Kratky-Porod concept |25j that accounts 
for interactions between nearest neighbors only, and its continuous analog, the model of worm- 
like chains [33| . 

Unlike most previous studies concerned with phantom chains, the present work concen- 
trates on constitutive equations for linear macromolecules modeled as self-avoiding random 
walks (SARW) |34j . Although an exact expression for the distribution function of end-to-end 
vectors of a SARW is unknown, two convenient approximations are widely used |35l I36j . Ac- 
cording to the first, the distribution function of end-to-end vectors is described by the stretched 
exponential expression 



po(Q) 



p° exp 



2,y 



(3) 



where Q = |Q|, and 0, 5, 
normalization condition 



R are positive constants. The pre-factor p° is found from the 
Po(Q)dQ = 1, (4) 



which implies that 



P 



1 

47T 



exp(-/?(|) 25 )Q 2 dQ 



The other approximation of the distribution function was developed in |371 138j 



Po(Q) 



exp 



(5) 



(6) 



In the literature on random walks |351 136j . Eq. (6) is conventionally referred to as the des 
Cloizeaux law. In the field of mathematical statistics, this formula is known [22] as a symmetrical 
Kotz-type distribution. Equation (6) involves four adjustable constants, a, /3, S and R. The 
pre-factor p° is determined from the normalization condition (4), 



R 



2n 



P 



4tt 



-p(-(!) 2 V i+ ^q 



(7) 



There is a substantial difference between the parameters a, (3 and R, on the one hand, and 
the exponent 5 in Eqs. (3) and (6), on the other. While the quantities a, (3 and R can adopt 
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practically arbitrary values, the parameter 5 is strongly connected [312 with the scaling exponent 
v in Eq. (2), 

5 = - (8) 
1 — v 

Due to some technical peculiarities of the renormalization-group method |40| employed for the 
evaluation of the scaling exponent, it is convenient to assess v for a SARW in a <i-dimensional 
space. Flory was the first who conjectured that 

d , . 

(9) 



d + 2 



For a three-dimensional space, Eqs. (8) and (9) imply that v = | and 

* 4 do) 

Although Eq. (9) has not been either proved or disproved up to now, the estimate (10) is used 
in the present study for two reasons: (i) formula (9) is widely employed in the field of statistical 
physics of polymers, and (ii) Eq. (10) allows necessary calculations to be performed explicitly. 
The objective of this study is three-fold: 

1. To derive an analytical expression for the strain energy W of a polymer chain with an 
arbitrary distribution function of end-to-end vectors Po(Q)- 

2. To apply this formula in order to develop explicit expressions for the strain energy of chains 
with distribution functions (3) and (6). 

3. To calculate the strain energy density of a network of self-avoiding chains, and to find 
adjustable parameters in the governing equations by fitting experimental data. 

To derive constitutive equations, we employ the Boltzmann formula, according to which the 
probability po(Q) that a chain has an end-to-end vector Q is expressed in terms of the configu- 
rational free energy t/o(Q) a s 

p (Q) = exp(-^M). (11) 

The function £7o(Q) determines the energy of a chain in the reference (stress- free) state. When 
macro-deformation is applied to a polymer network, the initial end-to-end vector Q of the chain 
is transformed into some vector q at an instant t > (time t = corresponds to the application 
of external loads). For an affine network of chains, the distribution function p(t, q) of the vector 
q obeys an appropriate Smoluchowski equation [HOI that is solved explicitly. 

Given p(t, q), the configurational free energy of a chain in the actual (deformed) state U(t, q) 
is described by the equation similar to Eq. (11), 



p(t,q)=exp(--^). 



:i2) 



When the functions t/o(Q) and U(t,q) are found, the increment of the configurational free 
energy is given by 

AU = U-U , (13) 
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The strain energy per chain W(t) is determined by averaging the increment AU with respect to 
an appropriate distribution function (j>q or p). A general expression for W will be developed in 
Section 3. At volume-preserving deformation of a Gaussian chain with the distribution function 



P0(Q) = (2^) 5eXP (-^) 



(14) 



where Q = |Q|, this formula is transformed into Eq. (1), when the averaging is performed 
with respect to the current distribution of end-to-end vectors, and it implies the Mooney-Rivlin 
constitutive law 



when the averaging is performed with respect to the functions p and po with the weights a and 
1 — a, respectively. 

Regrettably, the strain energy of a chain cannot be expressed as an elementary function of 
principal stretches for an arbitrary three-dimensional deformation. To develop such an ana- 
lytical expression, we apply an approximation procedure similar to the three-chain hypothesis 
|26| . according to which "the network of chains. . . is mathematically equivalent to three indepen- 
dent sets of chains respectively parallel to the axes of a three-dimensional Cartesian coordinate 
system" 

By using the three-chain approximation, explicit formulas are derived for the strain energy 
of self-avoiding chains with distribution functions (3) and (6). An important conclusion of our 
analysis is that for an incompressible network, the strain energy of a chain with the stretched 
exponential distribution function (3) is described (up to a small correction term) by the three- 
term Odgen law 



Unlike the original Ogden formula |43| . where a, a n and n n (n = 1, 2, 3) are treated as adjustable 
parameters, our expression contains only 2 material constants, while the exponents a n equal 1, 
3 and 5, respectively. This result provides a micro-mechanical basis for the Ogden model, on the 
one hand, and allows the number of experimental constants in Eq. (16) to be reduced noticeably, 
on the other. 

For the des Cloizeaux law (6), an expression for the strain energy per chain W is developed for 
an arbitrary (positive and negative) exponent a. The function W differs from that described by 
Eq. (16) by an additional term that depends on a. Surprisingly, approximation of experimental 
data reveals that the best fit of observations is reached at negative values of a. As these values 
have no physical meaning within the model of self-avoiding random walks (although they are 
not forbidden from the physical standpoint), we conclude that the Ogden constitutive law (16) 
describes the elastic response of all chains that can be modeled as self-avoiding random walks. 

Finally, the mechanical response of self-avoiding chains is analyzed at small deformations. 
At uniaxial tension (compression) with small strains, the elastic behavior of a chain is entirely 
described by the only parameter, an analog of the Young's modulus. The elastic modulus is 
proportional to k B T, but the coefficient of proportionality depends on the distribution function 
Po(Q)- It is found that self-avoiding chains with the stretched exponential distribution func- 
tion (3) are "stronger" than Gaussian chains with the distribution function (14) by a factor of 
three, whereas the ratio of the Young's modulus of a self-avoiding chain with the des Cloizeaux 
distribution function to that of a Gaussian chain linearly increases with exponent a. 



W = -feBT[a(Ji - 3) + (1 - a)(J 2 - 3) 



(15) 



1 r St -Ha i 

W = -k B T $> n £ [ a ( A "? - 3 ) + ( X - a X A ™ 2 - 3) . 



71=1 771=1 



(16) 
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The exposition is organized as follows. The Smoluchowski equation for the distribution 
function p(t, q) is resolved in Section 2 for an arbitrary time-dependent deformation of an affine 
network. A formula for the strain-energy of a chain is developed is Section 3. The strain 
energy of a Gaussian chain is calculated in Section 4. A three-chain approximation procedure 
is introduced in Section 5. The mechanical energy of a self- avoiding chain with the stretched- 
exponential distribution function is found in Section 6. This expression is applied in Section 
7 to fit experimental data on several elastomers. The strain energy of a chain whose statistics 
is governed by the des Cloizeaux law (6) is determined in Section 8, and the corresponding 
constitutive law is validated in Section 9. The elastic moduli of self-avoiding chains are calculated 
in Section 10. Some concluding remarks are formulated in Section 11. To avoid technical details 
in the main text, necessary calculations are given in Appendices. 



2 Transformation of the distribution function 



With begin with the solution of the Smoluchowski equation for the distribution function of end- 
to-end vectors p. Consider a chain with an end-to-end vector Q in the reference state and an 
end-to-end vector q is the actual state at time t > 0. In an affine network, transformation of 
the reference state of the chain into its deformed state is described by the formula 



q = F(t)-Q, 



(17) 



where F stands for the deformation gradient for macro-deformation (at this stage of the analysis, 
we do not impose the incompressibility condition on the tensor F). The function F(i) obeys the 
differential equation 

dF 

— = L-F, F(0) = I, (18) 

where L(t) is the velocity gradient, I is the unit tensor. 

The distribution function p(t, q) of the end-to-end vector q satisfies the equation [3U] 



dp 
dt 



d_ 
<9q 



"777 = -77- • L • qp 



-Xi(D)p-g-L.q, 



with the initial condition 



Here 



p(0,q) =p (q)- 



D = i[L + L' 



(19) 
(20) 



is the rate-of-strain tensor, T denotes transpose, and I m stands for the mth principal invariant 
of a tensor (m = 1, 2, 3). Simple algebra implies that the solution of Eqs. (19) and (20) reads 



p(t,q) =p (F- 1 (t)-q)exp[- j I^{p{s))ds 



(21) 



It follows from Eq. (18) that the third principal invariant of the deformation gradient obeys the 
equation 

dl 3 (F) 



dt 



Ti(D)J 3 (F). 



The solution of Eq. (22) with the initial condition X3(F(0)) = 1 reads 



X 3 (P(t)) = exp 



.(dm) 



ds 



(22) 



(23) 
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Equations (21) and (23) imply that 



P(t, q) = 



p (F- 1 (t)-q) 



(24) 



(25) 



urn) ' 

To demonstrate that function (24) satisfies the normalization condition 

J p(*,q)^q = i, 

we substitute expression (24) into the left-hand side of Eq. (25) and find that 

J Pfoq) d q = T3( F (t)) / Po(F- 1 (i)-q)cZq. 

Introducing the new variable Q by Eq. (17) and bearing in mind that 

dq = J 3 (F)dQ, 

we arrive at the formula 

J p(t, q)dq = J Po (Q)dQ. 
Equation (25) follows from this equality and Eq. (4). 



3 Strain energy of a chain 

Our aim now is to calculate the strain energy of a chain in an afiine network by using Eq. (24). 
The configurational free energies of a chain in the reference and actual states are connected with 
appropriate distribution functions by Eqs. (11) and (12), 



(26) 



U (Q) = -k B T\np (Q) 



U(t,q) = -fcnTlnpftq). 



(27) 



There are two ways to determine the strain energy per strand W. According to the first, the 
increment AU of the configurational free energy caused by transition from the reference state 
to the actual state, 

AU(t, Q) = U(t, Q) - U (Q) = -k B T[\np(t, Q) - lnpo(Q)' 
is averaged with the help of the distribution function in the reference state, 

W 1 (t) = -k B T J []np(t, Q) - ln Po (Q)]po(QMQ. (28) 

According to the other approach, we calculate the increment of the configurational free energy 
with respect to the actual state, 

Atf(i,q) = U (q) - U(t,q) = -fc B r[lnpo(q) - hxp(t, q) 
and average it by using the distribution function in the deformed state, 

W 2 {t) = -k B T f flnpo(q) - Iiip(t,q)lp(t,q)dq. (29) 
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We begin with transformation of Eq. (28). Substituting expression (24) into this equality and 
using Eq. (26), we obtain 



W^t) = -k B T\ / flnpotF- 1 ^) • Q) - lnpo(Q)WQ)dQ - InJ 3 (F(t)) 



(30) 



We now proceed with Eq. (29), combine this equality with Eq. (24), we find that 



W 2 (t) 



J [lnpotF- 1 ^) • q) - lnpo(q) - In l 3 (F(t))] Po (F-\t) • q)dq. 



Introducing the new variable Q by Eq. (17) and using Eqs. (25) and (26), we find that 



W 2 (t) = k B T{J [lnpo(Q) - lnpo(F(t) • Q)Jp (Q)dQ - lnJ 3 (F(t))|. (31) 

It seems natural to define the strain energy per strand W as the weighted sum of the strain 
energies W\ and W 2 that are determined by employing different ways of averaging of the con- 
figurational free energy, 

W = (l-a)W 1 + aW 2 , (32) 
where a 6 [0, 1] is a material parameter. Combining Eqs. (30)-(32), we arrive at the formula 

W{t) = A; B r|-(2a-l)lnJ 3 (F(t)) + y a(lnj>o(Q) - lnpo(F(t) ■ Q) 



+ 



(l-a)(lnp (Q)-lnpo(F- 1 (t)-Q)) 



Po(Q)dQ\. 



(33) 



To proceed with transformations of Eq. (33), an additional hypothesis is needed regarding the 
function po{Q)- Assuming the end-to-end vectors to be distributed uniformly in the reference 
state, we set 

Po(Q) = P(Q 2 ), (34) 

where Q 2 = Q • Q, and P(r) is a given function of a scalar argument r. Substituting expression 
(34) into Eq. (33) and taking into account that 



(F Q) (f Q) = Q F T F Q = Q C Q, 



= Q B 1 Q, 



2|(F) = Js, 



where the left and right Cauchy-Green deformation tensors are determined by the conventional 
formulas 



B = F • F , C = F • F, 



(35) 



we find that 



W 



A; B r|-^(2a-l)lnJ 3 + J a(lnP(Q 2 ) - lnP(Q • C • Q 



P(Q 2 )dQ 



+ (1 - a)(lnP(Q 2 - lnP(Q • B" 1 • Q)) 
It is easy to show (Appendix A) that 

J lnP(Q • B" 1 • Q)P(Q 2 )dQ = J lnP(Q • C" 1 • Q)P(Q 2 )dQ. 



(36) 



(37) 
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It follows from Eqs. (36) and (37) that the strain energy per strand reads 



W = fc B T|-i(2a-l)lnJ 3 + J 



+ (1 - a)(lnP(Q 2 ) - lnP(Q • CT 1 • Q) 



a(lnP(Q 2 )-lnP(Q-C-Q 



P(Q 2 )dQ . 



(38) 

Our aim now is to apply Eq. (38) in order to determine the strain energy of a Gaussian chain 
with distribution function (14). 

4 Strain energy of a Gaussian chain 

According to Eqs. (14) and (34), the function P(r) reads 

P (r) = (_^)* exp( _^). 

Simple algebra (Appendix B) implies that 

J (lnP(Q 2 ) - lnP(Q • C • Q))p(Q 2 )dQ = ±{J X - 3), 

J (lnP(Q 2 ) - mP(Q • CT 1 • Q))p(Q 2 )dQ = i(J_ x - 3), (39) 

where J_ m = X m (C _1 ) (m = 1, 2, 3). Substitution of expressions (39) into Eq. (38) implies that 



W = \k B T 



a(Ji - 3) + (1 - a)(J_i - 3) - (2a - 1) In J 3 



Bearing in mind that J_i = J2/J3, we present Eq. (40) in the form 



W 



1 



;k B T 



a{ Ji - 3) + (1 - a)(- 2 - - 3) - (2a - 1) In J 3 
^3 



(40) 



(41) 



Equation (41) implies that for an incompressible polymer network with J3 = 1, the strain energy 
of a chain is given by the Moony-Rivlin law (15). Equation (41) is reduced to the constitutive 
relation (1) for an incompressible neo-Hookean medium (which is traditionally associated with 
the response of Gaussian chains) at a = 1 only (when the averaging of the increment of config- 
urational free energy AU is performed with respect to the distribution function p in the actual 
state). It is worth noting that our results are in accord with those developed in [21 by using 
another approach (a tube model for rubbery polymers). 



5 Three-chain approximation 

Our aim now is to introduce a three-chain approximation for the strain energy per chain based 
on expression (40) for the mechanical energy of a Gaussian chain. Bearing in mind that 



J\ — Ai + A2 + A3 



J - — h— h — 

Ai A2 A3 ' 



J 3 — AiA2A 3 , 



where A m (m = 1,2,3) are eigenvalues of the Cauchy-Green tensor C, we find from Eq. (40) 
that 



W 



lk B T L(Ai + A 2 + A 3 - 3) + (1 - a) (— + — + — - 3) - (2a - 1) ln(Ai A 2 A3) 
I L \Ai A2 A3 ' 



(42) 
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Let us consider now three networks of chains. The principal axes of the Cauchy-Green deforma- 
tion tensor C^" 1 ) for the mth network (m = 1, 2, 3) coincide with the eigenvectors of the tensor 
C, whereas the corresponding eigenvalues {A^ , \% ■> A3 } are given by 



(1) 


— Ai, A^ 


= 1, 




= 1, 


(2) 


= 1, A< 2) 


= A 2 , 


X (2) 


= 1, 


(3) 
1 


= 1, Af 


= 1, 


X (3) 
A 3 


= A 3 



A 
A 
A 

It follows from Eqs. (42) and (43) that the strain energy per chain in the mth network reads 



(43) 



W {m) = K bT U Xm _ x) + (l _ a) - l) - (2a - 1) In A 



1 

Am 



Equations (42) and (44) imply that 



w= w(m) - 



(44) 



(45) 



m=l 



Formula (45) may be referred to as the three-chain approximation of the strain energy W (it 
slightly differs from the three-chain approach proposed in 26 ). Equation (45) provides an exact 
expression for the strain energy of a Gaussian chain. For non-Gaussian chains, this formula may 
serve as a convenient approximation for the strain energy. 

Returning to general expression (38) and applying the three-chain approximation, we find 
that the strain energy W is given by Eq. (45), where read 



A; B T|-i(2a-l)lnA m + J a(lnP(Q 2 ) - lnP(Q • C< m > • Q; 
+(1 - a)(lnP(Q 2 ) - lnP(Q • (C^)" 1 • Q))jp(Q 2 )dQ 



(46) 



We now choose a spherical coordinate frame {Q, 8, (ft}, whose axes coincide with the eigenvectors 
of the tensor cK m ) and the z-axis is directed along the eigenvector with the eigenvalue that differs 
from unity. The quadratic forms Q • ■ Q and Q • (C^)" 1 • Q are calculated as follows: 

Q • C<™> • Q = Q 2 (sin 2 9 + X m cos 2 0' " ~ ' ~ ' ' * ' ' 



Q-(C( m ))- 1 -Q = Q 2 (sin 2 ^ + 



A,, 



■ cos 



Substituting these expressions into Eq. (46) and performing integration over (ft, we find that 



k B T{ --(2a- l)lnA m + 2vr / P(Q 2 )Q 2 dQ sm6d9 



(in P(Q 2 ) - In P (Q 2 (sin 2 9 + A m cos 2 1 
+ (1 - a)(lnP(Q 2 ) - lnP(Q 2 (sin 2 9 + A" 1 cos 



Introducing the new variable z = cos 9, we obtain 



k B T 



1 



(2a - 1) In A m + 4vr / P(Q 2 )Q 2 dQ 



x jf 1 a(lnP(Q 2 ) - lnP(Q 2 (l + (A m - l)z 2 
+ (1 - a)(lnP(Q 2 ) - lnP(Q 2 (l + (A m x - l)z 2 ) 



dz 



(47) 



Our aim now is to calculate the strain energy W of a self-avoiding chain with distribution 
function (3) by using Eqs. (45) and (47). 
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6 The stretched exponential distribution function 

Equations Eqs. (3) and (34) imply that 



(48) 



where the pre-factor p° is given by Eq. (5). Substitution of Eq. (48) into Eqs. (45) and (47) 
results in (see Appendix C for detail) 



W 



k B T\--(2a-l)ln(X 1 X 2 X 3 ) + 



+ (1 " a) E + ^A- + - j- 6 + W) 

m=l 



a 

m=l 

5 

16 



16 



Am 



where 



E(«G(A m ) + (l-a)G(A- 1 ))-3 



m=l 



(49) 



(50) 



The function G(z) in Eq. (50) reads 

G(z) = < 



(1 — z) 2 arcsin(l — z) 2 , z < 1, 

1, ; , * =1 ' 

(z - 1)~2 ln(z2 + (z - 1)2), z > 1. 



(51) 



This function does not exceed § for any z > 0, and it monotonically decreases (however, rather 
weakly) with z and vanishes as z — ► 00. 

At finite strains, when the norm of the Cauchy-Green tensor C is large compared with unity, 
the last term in Eq. (49) may be disregarded. Under this assumption, Eq. (49) is reduced to 
the constitutive equation of the Ogden medium with a special choice of material parameters. 



W 



k B T 



1 ^"|5P\3^1"11 

-(2a - 1) In J 3 + a V (-Xl + — X m + — X m ) 

2 6 24 16 16/ 
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For an incompressible network, the first term in this equality vanishes, and we obtain, in accord 
withEq. (16), 
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(52) 



To demonstrate that Eq. (52) provides quite an acceptable approximation of the strain energy 
W for conventional loading programs, we calculate the strain energies determined by Eqs. (49) 
and (52) for uniaxial extension of an incompressible network, 



Ax = k 2 , A 2 = k~\ A 3 = k~\ 



(53) 
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where k stands for elongation ratio, and for simple shear with 




where k denotes shear. For definiteness, we set a = 0.5. The quantities W and Wq are plotted 
versus k in Figures 1 and 2. These figures show that the dimensionless strain energies W = 
W/(kB,T) practically coincide when they are calculated with and without the correction term 
W. 

Formula (52) determines the strain energy of an individual chain. Neglecting the energy of 
interaction between chains (this energy is conventionally accounted for by the incompressibility 
condition 30 ) , we calculate the strain energy per unit volume of a network as the sum of strain 
energies of chains, 

WM 5 A 5 ^l u \ 5-| 5-1 ll\ . A . 

w o = wL (e A - + + y^ a - - iq) + ^2 2. (e m + 24 m + le m ~ 16 ' ( } 

m=l m=l / 

where 

m = k B aMT, /j 2 = k B (l - a)MT, (55) 

and M is the number of chains per unit volume. Equation (54) involves only two coefficients, n\ 
and fi2, to be found by fitting observations. Our aim now is to show that constitutive law (54) 
correctly describes experimental data at uniaxial tension, uniaxial compression and equi-biaxial 
tension of elastomers. 



7 Results of numerical simulation 



At three-dimensional deformation of an incompressible medium without rotations, the principal 
Cauchy stresses a m are expressed in terms of the principal stretches A m by the formulas 



<7,< 



A, 



dw 



(56) 



where p is an unknown pressure, and w is a strain energy per unit volume. Substitution of 
expression (54) into Eq. (56) results in 



A r, 



/iiF (A 

r 



p- 



(57) 



where 



^o(A) 



5/1,3 1 1 1, 1. 

l(s A2 + 4 A2+ ^ 2) - 



(58) 



For uniaxial extension with an elongation ratio k, the principal stretches are given by Eq. (53). 
Excluding the pressure p from the condition &2 — &3 — 0, we find from Eq. (57) that the 
engineering stress a e = a\jk is given by 



Vi(k 2 F (k 2 ) - ifb(i)) + » 2 (kF (k) - ^F (^) 



1 



1 



(59) 



For equi-biaxial deformation of an incompressible material, the principal stretches A m read 



A, 



A 2 = k 



\ 3 = k~ 



(60) 
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where k stands for elongation ratio. We substitute Eq. (60) into Eq. (57), exclude the unknown 
pressure p from the condition 03 = 0, and find the Cauchy stresses a\ = 0% = a. The engineering 
stress <r e = a/k is determined as 



ae= k 



1 r /*i(*W) - +^(k i F (k^) - pfb(p))]. (6i) 



Given an experimental dependence a e (k), the coefficients /ii and //2 in Eqs. (59) and (61) are 
found by the least-squares technique from the condition of minimum of the function 



fen 

where the sum is calculated over all elongation ratios k n at which measurements are reported, 
<r| xp is the stress measured in a test, and <7g Um is given by Eqs. (59) and (61). 

We begin with fitting the observations on uniaxial extension of natural rubber with var- 
ious amounts of cross-linker tetramethylthiuram disulfide (TMTD) reported by Kliippel and 
Schramm The experimental data together with their approximation by Eqs. (58) and (59) 
are depicted in Figure 3 (where phr means parts per hundred parts of rubber). This figure 
shows that Eq. (59) with the only parameter [i\ (we set \ii = to reduce the number of mate- 
rial constants) provides good agreement between the observations and the results of numerical 
simulation. The elastic modulus [x\ is plotted versus concentration of cross-linker <f> in Figure 4, 
where the experimental data are approximated by the linear dependence 

Ml = tfV- (62) 

The coefficient /i^ in Eq. (62) is determined by the least-squares technique. Formula (62) 
confirms our assumption that the energy of interaction between chains in a network may be 
taken into account with the help of the incompressibility condition. It implies that the modulus 
/ii is proportional to the amount of cross-linker, which, in turn, is proportional to the number 
of chains per unit volume M, in accord with Eq. (55). 

We proceed with matching the experimental data on uniaxial compression of carbon black- 
filled chloroprene rubber reported by Bergstrom and Boyce pQ. The experimental dependencies 
a e (k) measured at various amounts of filler are depicted in Figure 5 together with their approx- 
imations by Eq. (59) with \ii = 0. The modulus [i\ is plotted versus concentration of filler in 
Figure 6. The experimental data for fJ>i((j>) are approximated by the linear dependence 

Mi=M? } +M?V, (63) 

where the coefficients /x^ an d Mi are found by the least-squares method. Figure 6 demonstrates 
that Eq. (63), which is conventionally employed in the mechanics of composites to describe the 
effect of filler on the shear modulus, correctly predicts the experimental data. 

Figures 3 and 5 show that the observations on natural and chloroprene rubbers at uniaxial 
tension and compression can be approximated by the model with a = 1, i.e. when averaging 
of the configurational free energy is performed by using the distribution function of end-to-end 
vectors in the deformed state. To demonstrate that this is not the rule in the general case, we 
fit experimental data on uniaxial tension-compression of carbon black-reinforced natural rubber 
reported by Hartmann et al. 9^ and of polybutadiene rubber presented by Roland et al. 
The observations and the results of numerical analysis are depicted in Figures 7 and 8, which 
reveal good agreement between the experimental data and the results of simulation. 
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To validate the governing equations, we approximate the observations on uniaxial extension 
of a synthetic rubber reported by Chevalier et al. find adjustable parameters fii and /i2, 
calculate the engineering stress a e at equi-biaxial tension by Eq. (61), and compare the results 
of numerical simulation with the experimental data. Figure 9 demonstrates that the model 
provides quite an acceptable prediction of the elastic response at equi-biaxial deformation with 
moderate finite strains, 1.0 < k < 1.4. 

Figures 3 to 9 reveal that Eq. (54) with two adjustable parameters can correctly describe the 
mechanical response of elastomers at large deformations. It is worth noting that this assertion is 
not trivial (the fact that the Ogden model adequately describes observations has been confirmed 
by numerous experimental data in the past three decades), because the number of material 
constants in Eq. (54) is noticeably smaller than that in the original expression (16) for the 
strain energy density. 

Another result to be mentioned is that the description of some experimental data by Eq. 
(54) is not perfect. As examples, we refer to the observations depicted in Figures 7 and 9. This 
conclusion may be attributed to the fact that Eq. (3) for the stretched-exponential distribution 
function of a self-avoiding random walk is overly simplified, and more sophisticated formula (6) 
should be used. Our aim now is to develop governing equations for a network of self-avoiding 
chains with the des Cloizeaux distribution function of end-to-end vectors. 



8 The des Cloizeaux distribution function 



For a SARW with the distribution function (6), the function P(r) reads 



exp 



r \S 



(64) 



where 5 is given by Eq. (10), and p° is determined by Eq. (7). Substitution of expression (64) 
into Eqs. (45) and (47) implies that (see Appendix D for detail) 



W = A; B r|-(2a-l)(a + i)ln(AiA 2 A 3 ) 

+ (l + r)[ a ^G Al+ 2l Al+ l| Ai 

m=l 
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-2a (o-H(Xm) + (1 - aMA™ 1 )) , 
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+ W 



where 



H{z) 
H(z) = 0, 
H{z) - ] 



arctan yjz— 1 



fz-1 



In 



l+Vl^z 



z > 1, 
z=l, 

1, Z<1. 



(65) 



(66) 



The function H(z) monotonically decreases in (0, oo) and tends to —1 when z — ► oo. The graph 
of this function is presented in Figure 10. 

Using notation (52) and neglecting the small term W compared with Wo, we find from Eq. 
(65) that 

3 



W=(l + ^a)W - k B T [(2a - l)(a + h In J 3 + 2a £ (aH(X m ) + (1 - a)H{\^) 



771=1 



(67) 
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In the tensor form, Eq. (67) reads (see Appendix E) 



W 



-(2a- l)(a + l)lnJ 3 + (l + |<* 
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where 



a(-l)%(E£) + (l-a)J 1 (Ei) 



Er^(C-I), E.\ 1(1 



2' " "2 
are the Cauchy and Almansi strain tensors, respectively. 

At volume-preserving deformations, Eqs. (52) and (67) imply that 



(68) 
(69) 



W 



= k B rL £ [(1 + |a) (1 \l + ^ A| + A a! - 11) - 2atf (A„ 

m=l 

+(1 - a) £ [(1 + ^)(^A m f + ^A m f + ^A m ^ - 11) - 2a#(A m 1 ) 
m i 



Applying the same approach that has been employed in Section 6, that is neglecting the energy 
of interaction between chains and assuming the strain energy of a network to coincide with the 
sum of strain energies of individual chains, we arrive at the formula for the strain energy density 
per unit volume of an incompressible network 



w = a 



m=l 

+(1 - a) £ [^i(^~ f + ^A~ f + _ 11) _ ^(A- 1 ) 



m=l 



where we set 



Mi = fcBT(l + -a)M, 



li 2 = 2k B TaM. 



(70) 



(71) 



In the general case, constitutive law (70) involves 3 adjustable parameters, a, fii and \i<i- It is 
worth noting a substantial difference between the parameter \i 2 in Eqs. (55) and (70). According 
to the physical meaning of the quantity \ii in Eq. (55), this parameter should always be positive. 
On the contrary, the parameter ^ is Eq. (70) may be positive [this corresponds to the case 
a > in Eq. (6)], as well as negative [when a 6 (— §, 0)]. The latter condition ensures that the 
integral of the function P(Q 2 ) over the entire space converges. 

Our aim now is to develop stress-strain relations for finite deformations of an incompress- 
ible network, to derive explicit expressions for the principal Cauchy stresses, and to compare 
experimental data on elastomers with results of numerical analysis. 



9 Fitting of observations 

It follows from Eq. (70) that 
dw r 

— - = a mF (X m ) + n 2 F(\ r , 
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HiF (—) + ii 2 F(—) 
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A, 



(72) 
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where Fq(z) is given by Eq. (58), and 

dH 



1 /arctanvz— 1 1 



/2-1 



F(z) 



dz 



1 

3' 

, 1 f 1 

2(1-2) V 2 2^ 



l n .i+v^ 



z > 1, 
3 = 1, 
2 < 1. 



(73) 



The function F(z) is positive and continuous at any z > 0, it monotonically decreases with z 
and vanishes when z — > oo. The graph of the function ^(-z) is depicted in Figure 10. 
Substitution of expression (72) into Eq. (57) results in the formula 



-p + HiX m aF (X 



1-a 



"A 2 " F ° l T 
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+ /i2A m aF(X 
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(74) 



According to Eq. (74), at uniaxial extension with an elongation ratio k, see Eq. (53), the 
engineering stress a e is given by 
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(75) 



At equi-biaxial extension with an elongation ratio k, see Eq. (60), the engineering stress reads 
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At biaxial extension with elongation ratios k\ and k 2 



\\ = kj, A 2 = k 2 , A 3 = (k\k 2 ) 2 , 
the engineering stresses cr e ,i = <T\jk\ and o" Cj 2 = o 2 jk 2 are determined by 
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(78) 



Our purpose now is to apply Eqs. (75), (76) and (78) in order to approximate observations at 
uniaxial, equi-biaxial and biaxial extension of elastomers. To find the constants a, /ii and fj, 2 , 
we divide the interval [0, 1], where the parameter a is located, into I subintervals by the points 
aW = i/\a (i = 1, — 1) with Aa = 1/L Given the coefficients fi\ and /i 2 are found 
by the least-squares method from the condition of minimum of the function R. The "best-fit" 
parameter a is determined from the condition of minimum of this function on the set {a^}. 

We begin with the experimental data on uniaxial tension-compression of two types of carbon 
black-filled natural rubber. The experimental stress-strain curves reported by Roland et al. [3^1 
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are depicted in Figure 11 together with their approximations by Eq. (75). This figure shows 
that Eq. (75) with a = 1 (averaging with respect to the distribution function in the actual 
state) and (12 < (a weakly-singular at Q = des Cloizeaux distribution function) ensures good 
agreement with the observations. 

The latter peculiarity of the model appear to be rather typical. To confirm this assertion, 
we approximate the experimental data on equi-biaxial extension of natural rubber reported by 
James at el. |46j . The observations together with their fit are plotted in Figure 12. It is worth 
noting some similarity between the material constants determined by matching the experimental 
data for three different types of natural rubber: the modulus fix is of order of (2 4) • 10~ 2 MPa, 
and the modulus \^\ is of order of 3 -f- 7 MPa. 

To validate the model, we analyze observations on biaxial extension of poly(dimethylsil- 
oxane) network (volume fraction of PDMS 4> = 0.463) reported by Kawamura et al. |47j . The 
adjustable parameters are found by matching the dependence 0^,2(^2) at k\ = 1.9. Afterwards, 
the functions cr e ,2(&2) at k% = 1.1, 1.3, 1.5 and 1.7 are predicted numerically, and the results 
of simulation are compared with the experimental data. Figure 13 shows excellent agreement 
between the observations and the model predictions. 



10 The elastic modulus of self-avoiding chains 

We return now to expressions (52) and (70) for the strain energy of an incompressible network 
of self-avoiding chains in order to assess the effect of excluded volume on the elastic modulus 
per chain at small deformations. 

For a network of Gaussian chains under uniaxial tension with an elongation ratio k, the 
engineering stress a e is given by the conventional formula 22 

cr e = (hi + /i 2 A; _1 )(fc - k~ 2 ), 

where 

Ml = -k B aTM, [i 2 = -k B (l - a)TM. (79) 

At small strains, when k = 1 + e with e < 1, we obtain a c = 3(/Ui + ^2)^- The Young's modulus 
of a network E is given by 

E=^. (80) 
e 

Substituting expressions (79) into Eq. (80) and introducing the Young's modulus per chain as 

E c = ^7, (81) 

M v ; 

we arrive at the formula 

E c = -k B T. (82) 

For a network with strain energy density (52), the engineering stress at uniaxial extension is 
given by Eqs. (58) and (59). Omitting simple algebra, we find that at small strains, 

305, 

°e = — (^1 +/i 2 )e- 
It follows from this equality and Eqs. (55), (80), and (81) that 

E c = ™k B T. 
64 
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Combining Eq. (82) and this equation, we arrive at the formula 



^ = 595-8.1771, (83) 

£Gauss 96 ^ J 

which means that the "rigidity" of a self-avoiding chain exceeds that of a Gaussian chain by 
a factor of three (approximately). A reason for the anesthetic ratio in Eq. (83) is that the 
contribution of W into the strain energy density W is neglected. 

We proceed with the evaluation of the Young's modulus for a chain whose distribution 
function is described by des Cloizeaux law (6). It follows from Eq. (73) that 

dz v ' 5 

Using this expression, we find from Eqs. (58) and (75) that at uniaxial extension with small 
strains, 

„ 305 1 

Substitution of expressions (71) into this formula results in 

/305 1333 \, _ 

E c = [ — + a)k B T. 

\ 64 480 / 

It follows from this equality and Eqs. (82) and (83) that 

pdcsCloizcaux pOgdcn i ooo 

= + -P-a « 3.1771 + 1.8514a. (84) 

^Gauss ^"Gauss 72O 

Equations (83) and (84) imply that the Young's modulus per chain E c is independent of the 
characteristic length R and the exponent (5 in Eqs. (3) and (6), and it increases with a in Eq. 
(6) being proportional to 2a (approximately). 

Qualitatively, the latter conclusion is not surprising. The growth of a is tantamount to the 
fact that the tail of distribution function (6) becomes thicker, while the latter implies (within the 
conventional "force-stretch" approach in the statistical mechanics of polymers) an enhancement 
of the rigidity of chains. 

To assess the Young's modulus found in our approximation of the observations, we use Eq. 
(71), which implies that 

1/1 Lt!\-l 

a = -2v3-^ • 

Using the results presented in Figures 11 to 13, we find that a = —1.46 (Figure 11) and —1.48 
(Figure 12) for natural rubber, and —1.49 (Figure 13) for PDMS. For these parameters, Eq. 
(84) implies that the ratio ^dcsCioizcaux/^Gauss ig doge to Q M for PDM s and is located between 
0.47 and 0.49 for natural rubber. The difference between these quantities and the ratio of elastic 
moduli predicted by Eq. (83) may provide an explanation, why constitutive law (52) does not 
provide an adequate approximation of the observations depicted in Figures 11 to 13. 



11 Concluding remarks 

A general expression has been derived for the strain energy of a polymer chain with an arbitrary 
distribution function of end-to-end vectors. This formula has been applied to determine the 
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strain energy of a chain that is modeled as a self-avoiding random walk. By using three-chain 
approximation method, we has found that the strain energy of a network of chains with the 
stretched exponential distribution function coincides with the Ogden law with a special choice 
of adjustable parameters. In this case, the stress-strain relations involve only two material con- 
stants that are found by matching experimental data at uniaxial tension, uniaxial compression 
and equi-biaxial tension of elastomers. Good agreement is demonstrated between the observa- 
tions and the results of numerical simulation. 

For the generalized stretched-exponential distribution function (the des Cloizeaux law), a new 
expression is developed for the strain energy density of a network. In this case, the governing 
equations contain three material constants to be found by matching observations. It is shown 
that the stress-strain relations provide quite an acceptable fit of observations at uniaxial and 
bi-axial extension of filled elastomers even when the Ogden relations poorly approximate the 
experimental data. It is demonstrated that the latter case corresponds to rubber-like materials, 
whose elastic modulus at uniaxial extension is lower than that for Gaussian chains. 
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Appendix A 

To derive Eq. (37), we, first, choose a spherical coordinate frame {Q, 0, </>}, whose axes coincide 
with the eigenvectors of the symmetric tensor B _1 . In the spherical coordinates, the quadratic 
form Q • B 1 • Q is given by 

Q B 1 Q = Q 2 [(.B _1 )i sin 2 cos 2 + (B -1 ) 2 sin 2 sin 2 + (B~ 1 ) 3 cos 2 

where {B~ l ) m (m = 1,2,3) are eigenvalues of the tensor B -1 . The integral on the left-hand 
side of Eq. (37) reads 

//•oo t"K r2n r 

lnP(Q • B" 1 • Q)P(Q 2 )dQ = / P(Q 2 )Q 2 dQ sin0d0 (5 _1 )i sin 2 cos 2 
Jo Jo Jo L 

+ { B ~ l h sin 2 9 sin 2 + (5 _1 ) 3 cos 2 #] 

We now choose the spherical coordinate frame {Q' ,0' whose axes coincide with the eigen- 
vectors of the symmetric tensor C -1 . Repeating the above calculations, we find the integral on 
the right-hand side of Eq. (37), 



/roo rn r2n r 

lnP(Q-C- 1 -Q)P(Q 2 )dQ = / P(Q' 2 )Q' 2 dQ' sm0'd0' (CT^i sin 2 0' cos 
JO JO JO L 

+ (C~ 1 ) 2 sin 2 0' sin 2 <p' + (C _1 ) 3 cos 2 0'] d<f>' , 

where (C~ 1 ) m (m = 1,2,3) are eigenvalues of the tensor C _1 . Bearing in mind that the 
eigenvalues of the Cauchy-Green tensors B _1 and C^ 1 coincide, we conclude from these relations 
that Eq. (37) is satisfied. 
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Appendix B 

We begin with transformation of the term 
J lnP(Q ■ C • Q)P(Q 2 )dQ = | In 



2vr6 2 26 2 V2tt6 2 



|(Q.C-Q)exp(-^)dQ. (B-l) 



In spherical coordinates {Q,0,4>} directed along the eigenvectors of the symmetrical tensor C, 
the integral in Eq. (B-l) is presented as 

y(Q-C-Q)exp(-^)dQ = J™ Q A ^ V (-^)dQ J\m0d0 J^[\ lS m 2 cos 2 (f) 

+A2 sin 2 sin 2 cf> + A3 cos 2 

where X m (m = 1, 2, 3) are eigenvalues of the tensor C. The integrals on the right-hand side of 
this equality are calculated explicitly, 



£V«p(-f£)*j-s(:)*( 

f W sm0d0 f 
Jo Jo 



3 



27T r T A-TT 

Ai sin cos </> + A2 sin sin + A3 cos d<j) = -^-Ji- 
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Substitution of these expressions into Eq. (B-l) results in 



J lnP(Q ■ C ■ Q)P(Q 2 )dQ = | In 



2vr6 2 



Bearing in mind that P(Q 2 ) = P(Q • I ■ Q), we find from this equality that 

J lnP(Q 2 )P(Q 2 )dQ = hn 
Equation (39) follows from Eqs. (B-2) and (B-3). 
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Appendix C 

Substitution of Eq. (48) into Eq. (47) implies that 

= k B T^-^(2a-l)\n\ m + K £ (a(l + (\ m - l)z 
+ (l-a)(l + (A- 1 -l)z 2 ) 5 )dz-l]}, 



2\<5 



(C-l) 



where 



K = 4ttPP exp -0 

Jo 1 K R 



R 



Combining this equality with Eq. (5) and introducing the new variable r = Q / R, we find that 

K _ p Sfexp(-Pr 3 )r 5 dr ^ 
Jq 00 exp(— [3r 3 )r 2 dr 

We set z = (33 r in this equality and obtain 

K 



Jq 00 exp(— z s )z 5 dz 



J °° exp(— z 3 )z 2 dz ' 
The ratio of the integrals is easily calculated by integration by parts, 

K = 1. 

Substitution of expression (C-2) into Eq. (C-l) results in 

= fc B Tj-i(2a-l)lnA m + £ (a(l + (A m - l)z 2 ) 5 
+(l-a)(l + (A m 1 -l)z 2 ) 5 ) ( iz-l]}. 



(C-2) 



(C-3) 



To find an analytical expression for the function W^ m \ it is necessary to calculate the integral 

L(ji,S) = [\l + nz 2 ) 5 dz (C-4) 
Jo 

for fi = X m — 1 and fi = A" 1 — 1, and for 5 given by Eq. (10). We begin with the case fi > 0. 
Integration by parts implies that 



/< 



1 + [iz 2 )2dz = i z(l + fiz 2 



) f +5/(1 



+ fiz 2 )2dz 
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It follows from this formula that 



(C-5) 



where the function L is given by Eq. (C-4). Applying integration by parts to the new integral, 
we obtain 

J (1 + fiz 2 )idz = - z(l + Liz 2 )i +3 j {1 + fiz 2 )^dz , 

which results in 



^>2 ) = 4 



(1 + A*)3 + 3L(m, 



It is easy to check that 

y (i +^ 2 )* 



2(fe = - 

2 



z^Jl + fiz 2 + — ln(^/Lz + ^1 + ^z 2 ) 
yM 



It follows from this equality that 



L(p, -) = - [VTTJ1 + -j= ln(VJZ + VI + A*) 
Combination of Eqs. (C-5) to (C-7) yields 

Hp, \) = ^(i + m) 1 + ^(i + m) 1 + ^(i + m)^ + ^= MVm + 

Returning to the initial notation, we find from Eqs. (C-4) and (C-8) that 
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(1 + (A m - \)z 2 )idz = \\i + ^\i + ^-\i + ^-(\ m -i)- 

b 24 lb lb 

xln(Al + (A m -l)5) (A m >l). 
We now consider the case jj, < 0, set Mi = — //, and calculate the integral 
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mz 2 )2dz (m > o). 



Repeating the above transformations, we find that 



Mmi, -) 



r 3 
(l-/xi)f +5Li(/ii,-) 

(1-Mi) 1 +3L 1 (/i 1 ,i) 

.i 1 - 
(1 — /ii) 2 -| — arcsin y^Tf 



To develop the latter equality, we use the formula 

1 



Mi£ ) 2 = - 1 — [i\z 2 H arcsin y^Mi^ 



Mi 



It follows from Eq. (C-10) that 



r / 5, 1 .., .5 5 ,., .3 5 ,., ,i 5 
Mmi, = g(l - Mi)" + ^(1 - + ^(1 - Mi)* + ^= 



arcsin y'/ii 



(C-9) 



(C-10) 



(OH) 
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Bearing in mind that \i\ = 1 — A m , we find from Eq. (C-ll) that 

f\l + (\m-±)z 2 )ldz = + + ^-Xl + A (l-A m )-5 

Jo fc> 24 lb lb 

x arcsin(l - A m )^ (A m < 1). (C-12) 

It is obvious that 

J\i + (\ m -i)z 2 )idz = \^ + YA x ^ + h^ n+ h ( A ™ = 1 )- ( c - 13 ) 

Substituting expressions (C-9), (C-12) and (C-13) into Eq. (C-3) and using Eq. (45), we arrive 
at Eqs. (49) to (51). 



Appendix D 



It follows from Eq. (64) that for any m = 1,2,3, 
lnP(Q 2 )-lnP(Q 2 (l + (A m -l)z 2 ) = 



Kff [(l + (A m -D,f-l 



-aln(l + (A m - l)z 2 ), 



lnP(Q 2 )-lnP(Q 2 (l + (A- 1 -l)z 2 ) = /?(|) 25 [(l + (A" 1 - ) -1 

— aln(l + (A" 1 — l)z 2 ). 
Substitution of these expressions into Eq. (47) implies that 

W (m) = ^(m) _ ^(m) 

with 



(m) _ 



fe B T<j--(2a-l)lnA m + K 



jf (a(l + (A m - l)z 
+ (l-a)(l + (A m 1 -l)z 2 )^dz-l]}, 



(m) 



kftTKi / a ln(l + (A m - l)z z ) + (1 - a) ln(l + (A~ - l)z 2 ) dz, (D-2) 
jo L J 



where we set 



K 



47T/3 



(D-l) 



/ P(Q 2 )Q 2 ( 1+<5 )(fQ, = 4^a / P{Q 2 )Q 2 dQ. (D-3) 

jo jo 



Introducing the new variable z = Q/R in the first equality in Eq. (D-3) and using Eq. (7) to 
evaluate the other equality, we find that 



K = 



/ °° exp(-pz 3 )z 2a+5 dz 
J oo exp(-^ 3 )z 2Q + 2 dz' 



K\ = a. 



(D-4) 



To determine the coefficient K, we set y = f3sz, which results in 



/ oo exp(-/3y3)y2a + 2 d2/ - 
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Calculating the integral in the nominator by parts, we arrive at the formula 



K = 1 + -a. 



(D-5) 



The function $[ m - ) has already been found in Appendix C, 



(m) _ 



^bT<{ --(2a - l)lnA m + A" 



a(^Al + ^Al + ^A| t )+(l-a) 



1-5 5-3 
\ 2 _i_ \ 2 
— Am ~T~ /\ Til 

6 24 

(m) 



+ i^) + ^( aG ( A ™) + ( J - °) G ( A -)) " X ] }• ( D - 6 ) 



To determine the function ^J, , we introduce the notation 



M(ji) = / 1 ln(l + fiz 2 )dz 
Jo 



(D-7) 



and calculate the function M(/x) for ^ = A m — 1 and A m x — 1, respectively. Performing integration 
by parts in Eq. (D-7), we obtain 



r -i2=l /■! dz r /"^ 

M( / u) = zln(l + ^ 2 )-2 +2/ — — =ln(l + M )+2 / t 
L J 2=0 Jo 1 + L J0 1 



(D-8) 



Calculating the integral, we find that 



This implies that 



M(fi) = ln(l + n) + 2 ( arct ^ - l) , /U >0, 
M(/x) = 0, M = 0, 

M( M )=ln(l + ^)+2(^lni±^-l), ^ < 0. 

/* ln(l + (A m - l)z 2 )dz = In X m + 2H(X m ), 
Jo 



(D-9) 



where the function -ff(z) is given by Eq. (66). Substitution of expression (D-9) into the second 
equality in Eq. (D-2) results in 



(m) 



ksTK^ [(2a - 1) In A m + 2(aH{\ m ) + (1 - a)H{X^ 



(D-10) 



It follows from Eqs. (D-l), (D-4), (D-5), (D-6) and (D-10) that 



k B T{ l)lnA T 





2 \ 


— a 




3 J 



/l 5 5 a 5 i\ 
a(-A^ + -A^ + -A^) + (l-a) 



(^ Am§ + 2l Am§ + l! Am ") + ^( aG ( A -) + " a ) G ( A ™)) " X ] } 



-fc B Ta [(2a - 1) In A m + 2(off(A m ) + (1 - a)F(A m 1 
Equation (65) follows from Eqs. (45) and (D-ll). 



(D-ll) 
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Appendix E 



According to Eq. (52), the function Wq is given by 



Wo 



k B T 



<^ l(Cl) + 2l Jl(Cl)+ ll5 Jl(C " 



•6 v ' 24 v ' 16 

+ (i-a)(I Il (c-l) + Ax 1 (c-i) + A Il(c -i))-| 



(EM) 



Expanding the function H(z) in Eq. (66) into a Taylor series in the vicinity of the point z = 1, 
we find that 



2fc + 1 



k=i 1 fc=i 

which means that for any z belonging to the region where the series converges, 

{l-z) k 



H{z) = J2 



k=l 



2k + 1 



(E-2) 



It follows from Eq. (E-2) that 

3 



E = E am Ed - a-)' = E - o'>- 

m=l fc=l m=l fc=l 

Similarly, 

3 oo i 

E^ m 1 ) = E^- T ^(( I - c - 1 ) fe )- 

m=l fc=l 

Substitution of these expressions into Eq. (67) results in 

2 1 
(l + -ajW - k B T[(2a - l)(a + -) In J 3 

OO -I 

+ 2 « E ^r( aJl((I " c)k) + (1 " a)Jl((I " c ~ 1)fc) ) 



W 



(E-3) 



fc=i 



Equation (68) follows from Eqs. (E-l) and Eq. (E-3) and definition (69) of the strain tensors. 

As the series in Eq. (E-3) converges rather slowly, it is tempting to truncate it up to terms 
of an order ko with respect to the norms of the strain tensors. Taking into account one term 
only (ko = 1), we find that 



k 

E 

k=i 



i 



2k + 1 



aX!((I - C) fc ) + (1 - a)Ji((I - C" 1 )*)) 



- 1 - [a(Xi(C) - 3) + (1 - a^TitC- 1 ) - 3) 



(E-4) 



It follows from Eq. (E-4) that in the first-order approximation, the account for the pre- 
exponential term in the distribution function (6) is equivalent to the introduction of additional 
Gaussian chains into the network. 
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Taking into account two terms (ko = 2) in Eq. (E-3), we find, after simple algebra, that 

1 



£ 

k=i 



2k + 1 



oTxKl - C)*) + (1 - o)Ji((I - C- 1 )^) 



= Tg{ a [ 3 ( J i " 3 ) 2 + 7 ( J i " 3 ) " 6 ( J 2 " 3 ) 

+ (1 - a) [3(J_i - 3) 2 + 7(J_i - 3) - 6(J_ 2 - 3)] |. 



(E-5) 



According to Eq. (E-5), in the second-order approximation, the account for the pre-exponential 
term in Eq. (6) results in a decrease in the strain energy of a self-avoiding chain. 
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prene rubber reinforced with carbon black (CB). Symbols: experimental data pQ. Unfilled 
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numerical simulation. 
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Figure 10: The functions H(z) (thin curve) and F(z) (thick curve). 

Figure 11: The engineering stress a e versus elongation ratio k at uniaxial tension-compression 
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natural rubber. Circles: experimental data |46| . Solid line: results of numerical simulation with 
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Figure 13: The engineering stress o e ^. versus elongation ratio k 2 at biaxial tension of PDMS 
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